arXiv:1508.02357v3 [gr-qc] 14 Sep 2016 


The PyCBC search for gravitational waves from 
compact binary coalescence 


Samantha A. Usman^’^, Alexander H. Nitz^’^, Ian W. Harry^’^, 
Christopher M. Biwer^, Duncan A. Brown^, Miriam Cabero^, 
Collin D. Capano^’®, Tito Dal Canton^, Thomas Dent^, 

Stephen Fairhurst^, Marcel S. KehF’^, Drew KeppeF, 

Badri Krishnan^, Amber Lenon^, Andrew Lundgren^, 

Alex B. Nielsen^, Larne P. Pekowsky^ Harald P. Pfeiffer^® ®, 
Peter R. Saulson^, Matthew West^, Joshua L. Willis^’^. 

^ Department of Physics, Syracuse University, Syracuse, NY 13244, USA 
^ Cardiff University, Cardiff CF24 3AA, United Kingdom 

^ Max-Planck-Institut fiir Gravitationsphysik, Albert-Einstein-Institut, D-30167 
Hannover, Germany 

Max-Planck-Institut fiir Gravitationsphysik, Albert-Einstein-Institut, Am 
Muhlenberg 1, D-14476 Golm, Germany 

^ Maryland Center for Fundamental Physics & Joint Space-Science Institute, 
Department of Physics, University of Maryland, College Park, MD 20742, USA 
® Canadian Institute for Theoretical Astrophysics, University of Toronto, Toronto, 
Ontario MSS 3H8, Canada 

^ Max-Planck-Institut fiir Radioastronomie, Auf dem Hiigel 69, 53121 Bonn, 

Germany 

® Canadian Institute for Advanced Research, 180 Dundas St. West, Toronto, ON 
M5G 1Z8, Canada 

® Abilene Christian University, Abilene, TX 79601, USA 
E-mail: SEmiantha.usman@ligo.org 

Abstract. We describe the PyCBC search for gravitational waves from compact- 
object binary coalescences in advanced gravitational-wave detector data. The search 
was used in the first Advanced LIGO observing run and unambiguously identified 
two black hole binary mergers, GW150914 and GW151226. At its core, the PyCBC 
search performs a matched-filter search for binary merger signals using a bank of 
gravitational-wave template waveforms. We provide a complete description of the 
search pipeline including the steps used to mitigate the effects of noise transients in 
the data, identify candidate events and measure their statistical significance. The 
analysis is able to measure false-alarm rates as low as one per million years, required 
for confident detection of signals. Using data from initial LIGO’s sixth science run, we 
show that the new analysis reduces the background noise in the search, giving a 30% 
increase in sensitive volume for binary neutron star systems over previous searches. 
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The detection of the binary black hole mergers GW150914 and GW151226 by the 
Laser Interferometer Gravitational-wave Observatory (LIGO) has established the held 
of gravitational-wave astronomy [H [21 E]- Advanced LIGO [H E] will be joined by the 
Virgo and KAGRA [8] detectors in the near future, forming an international 

network of gravitational-wave observatories. Beyond the expected regular detections of 
binary black holes [9], compact-object binaries containing a neutron star and a black 
hole or two neutron stars are likely candidates for detection by this network [ininiE2]. 
Gollectively, these three types of gravitational-wave sources are referred to as compact 
binary coalescences (GBG). 

The identihcation of gravitational-wave candidates in the detector data is a complex 
process, which is performed by a search pipeline. The search pipeline is also responsible 
for determining the signihcance of each identihed gravitational-wave event. This paper 
describes a new pipeline to search for gravitational waves from compact-object binaries 
that uses the PyCBC software framework [13] to implement an all-sky search for 
compact binary coalescence. The PyGBG search pipeline described here builds upon the 
algorithms miiisiiiniiniiiEi used to search for compact-object binary coalescence with 
the hrst-generation LIGO and Virgo detectors jlHl [20l [211 122| [231 IMl [SSI [26l [2I|. This 
pipeline incorporates new algorithms that improve the sensitivity of the search and that 
reduce its computational cost. We provide a complete description of the search pipeline, 
emphasizing the new developments that have been made for the advanced-detector era. 
To demonstrate the efficacy of PyGBG, we re-analyze data from Initial LIGO’s sixth 
science run and show that the new pipeline can achieve a ~ 30% increase in sensitive 
volume to binary neutron stars, as compared to the pipeline used in the original analysis 
of the data [27|. Details of the results of the PyGBG search on the hrst aLIGO observing 
run are given in [2E1 [12] with the two observed binary black hole mergers discussed in 
detail in HIE]. 

This paper is organized as follows: Section provides an overview of the search 
pipeline and the methods used to detect gravitational waves from compact-object 
binaries in LIGO data. Section gives a description of the developments implemented 
in this pipeline. Section compares the performance of the new pipeline to that of the 
pipeline that analyzed the sixth LIGO science run and Virgo’s second and third science 
runs m Finally, Section summarizes our hndings and suggests directions for future 
improvements. 

2. Search Overview 

The purpose of the PyGBG search is to identify candidate gravitational-wave signals 
from binary coalescences in the detector data and to provide a measure of their statistical 
signihcance. Since the amplitude of the majority of gravitational-wave sources will 
be comparable to the noise background, signal processing techniques are required to 
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identify the candidate events. The gravitational waveforms for compact-object binaries 
can be modeled using a combination of analytical and numerical methods |29l l30] . 
Consequently, it is natural to use matched filtering to distinguish signals from the 
noise background [3T|. If the detector output contained only a stationary, Gaussian 
noise, the matched filter signal-to-noise ratio (SNR) would suffice as a detection statistic 
since its distribution in stationary, Gaussian noise is well known |32]. In practice, the 
detector data contains non-stationary noise and non-Gaussian noise transients |331 El] ■ 
Additional steps must therefore be taken to mitigate the effect of these noise transients 
and to assign an accurate statistical significance to candidate signals. Even for loud 
sources, the statistical significance of candidate detections must be empirically measured 
since it is not possible to shield the detectors from gravitational-wave sources and no 
complete theoretical model of the detector noise exists [S3] . 

Figure [T] shows the steps that the pipeline uses to find signals and measure their 
significance. The input to the pipeline is the detector’s calibrated strain data [36l 137] . 
In addition to possible signals, the strain data contains two classes of noise: a primarily 
stationary, Gaussian noise component from fundamental processes such as thermal 
noise, quantum noise, and seismic noise coupling into the detector [3] and non-Gaussian 
noise transients of instrumental and environmental origin [38l IM] . To eliminate the 
worst periods of detector performance, data quality investigations [331 SHI IHl ES] are 
used to characterize detector data into three general classes: (i) the data is polluted 
with enough noise that the data should be discarded without being searched, (ii) 
the data can be filtered, but candidate events that lie in intervals of poor data 
quality should be discarded, or vetoed, due to the presence of a instrumental or 
environmental artifacts, or (iii) the data is suitable for astrophysical searches. Data 
quality investigations are conducted independently of the search pipeline (by looking 
only at detector performance), as well as by determining the effect of instrumental 
artifacts on the noise background of the search pipeline. After removing data that is 
not suitable for astrophysical searches, the pipeline begins its analysis, following the 
remaining steps shown in Figure 

We do not a priori know the parameters of gravitational waves in the data, 
so a bank of template waveforms is constructed that spans the astrophysical signal 
space jia ESI 0111131 ESI igilHl SHI EHlEIlEaEs]. if the total mass of a compact-object 
binary is lower than M < 12 Mq [SH 03] and the dimensionless angular momenta of 
the compact objects (their spin) are small cSi^ 2 /Gm\ 2 0.4 [3Sl [37] (as is the case 

for binary neutron stars, where mi and m 2 and Si and S 2 are the component masses 
and spins respectively), then the detectors are sensitive to the inspiral phase of the 
waveform and this can be well modeled using the post-Newtonian approximations (see 
e.g. Ref. [23] for a review). For high-mass and high-spin binaries, analytic models 
tuned to numerical relativity can provide accurate predictions for gravitational waves 
from compact binaries |3Hl EHl [ 6 OI EOl [61] . The template bank is constructed to cover 
the space of circular binaries with aligned spins. It is generated so that the loss in 
matched-filter SNR due to the discrete nature of the bank is no more than 3%. The 
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Figure 1. A flowchart indicating the different steps of the search pipeline. Data 
from the detectors are averaged to create a power spectral density necessary to place 
a bank of templates that cover the search parameter space (blue box). Times when 
the detector data contain loud noise transients are removed, or vetoed. The data from 
each detector is then matched filtered and triggers are generated by thresholding and 
clustering the signal-to-noise ratio time series. A chi-squared test is computed for each 
trigger and the trigger’s matched-filter SNR is re-weighted by the value of the chi- 
squared statistic to better distinguish between signal and noise (yellow boxes). The 
pipeline determines which triggers survive time and template coincidence, discarding 
triggers that lie in times of poor data quality (red box). The triggers that pass the 
coincidence and data quality tests are labelled candidate events. Finally, multiple 
time shifts help generate a noise background realization that is used to measure the 
significance of the candidate events (bottom box). 
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exact placement of the templates depends on the detector’s noise power spectral density 
(PSD) Sn{f)- The PyCBC search pipeline places templates using a single noise PSD, 
averaged over all of the time and all detectors in the search network. The data from 
each detector in the network is matched hltered against this single template bank. We 
describe the process of constructing this template bank in Section 3.1[ 

Since the waveform of the target signals are well modeled, the pipeline uses matched 
hltering to search for these signals in detector noise. In the PyCBC pipeline, each 
detector’s data is hltered independently. The search templates are restricted to spin 
aligned binaries, and only the dominant gravitational wave harmonic [62]. Consequently, 
the sky location and orientation of the binary affect only the overall amplitude and phase 
of the waveform. These are maximized over when constructing the matched hlter by 
projecting the data signal against two orthogonal phases of the template h{t), given by 
hcos and hgin HT]. The matched hlter then consists of a weighted inner product in the 
frequency domain used to construct the signal-to-noise ratio (SNR), p(t), as: 


p\t) = 


(s I hcos) 

(hcos I hcof 


+ 


(s|hs 


(hsin I hsin) 


(s I hcos) T ("5 1 hsin) 

(hcos I hcos) 


( 1 ) 


where the inner product is given by 
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Here s(/) denotes the Fourier-transformed detector data, dehned by 

^-1-00 

s(/) = / 


( 2 ) 


(3) 


and h(/) denotes the Fourier-transformed template waveform. Sn{f) is the one-sided 
PSD of the detector noise dehned by 


(5(/)i(/')) = lsM)Hf - /'), 


(4) 


where angle brackets denote averaging over noise realizations and 6 is the Dirac delta 
function. The frequency limits /low and /high are determined by the bandwidth of the 
detector’s data, and the two phases of the template are related by hsin(/) = ihcosif)- 
Despite extensive data quality investigations, noise transients of unknown origin 
still remain after data quality vetoes are applied to the search. To mitigate the ehect of 
these noise transients, the pipeline identihes excursions in the input strain data s{t) 
and then applies a window to the detector data, zeroing out the data around the 
time of a noise transient before hltering. This procedure, called gating, is described in 


Section 3.2 Having removed these noise transients, the pipeline computes the matched- 
hlter SNR p{t) for each template in each detector. The search identihes the times when 
the matched-hlter SNR exceeds a predetermined threshold for a given template in a 
detector’s data. The pipeline applies a clustering algorithm, which takes the largest 
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value within a predefined window of p{t) and identifies maxima in the SNR time series. 
This process yields a list of times when a signal may be present, which are called 
triggers. The matched filtering, thresholding and clustering algorithms are described in 
Section 13.31 

Triggers generated by matched filtering the data against the template bank are 
subject to a chi-squared test that determines if the time-frequency distribution of 
power in the data is consistent with the expected power in the matching template 
waveform [16]. To construct this test, the template is split into p frequency bins. These 
bins are constructed so that each contributes an equal amount of power to the total 
matched-filter SNR. The matched-filter SNR, pi, is constructed for each of these p bins. 
For a real signal, pi should contain 1/p of the total power. The statistic compares 
the expected to the measured power in each bin according to 


=V^ 



(5) 


where p/.^^ and are the SNRs of the two orthogonal phases of the matched filter. 
Lower-mass binary systems, such as binary neutron stars, lose energy to gravitational 
waves more slowly than higher-mass systems. Consequently, the waveforms of lower 
mass systems are longer, having more gravitational-wave cycles in the sensitive band 
of the detector. The PyCBC pipeline allows the number of bins to be specified as a 
function of the intrinsic parameters of the template. This allows the search to use more 
bins in the chi-squared test for longer templates, making the test more effective. In 
previous analyses, the chi-squared test was the most computationally costly part of 
the pipeline [18]. The PyCBC pipeline uses a more efficient algorithm for computing 
the chi-squared statistic, which vastly reduces the computational cost. Details of the 
chi-squared test are given in Section 3.4] 

For a trigger of given matched-filter SNR, larger values of indicate a higher 
likelihood of a noise transient origin as opposed to a signal. For signals, the reduced 
chi-squared, y^ = y^/(2p —2), should be near unity. To down-weight triggers caused by 
noise transients, the matched-filter SNR is re-weighted [2ZIIIH] according to 


P = 


p/[(l + (Xr)^)/2]® , ifXr>l, 
p, if Xr < 1- 


( 6 ) 


Having computed the re-weighted SNR for each trigger, the pipeline discards all triggers 
that lie below a pre-determined re-weighted SNR threshold. 

The search requires that signals are observed with consistent parameters in the 
detector network. First, any triggers that occur during times of instrumental or 
environmental artifacts, as determined by the input data quality metadata, are vetoed. 
To be considered a candidate event, triggers must be observed with a time of arrival 
difference less than or equal to the gravitational-wave travel time between detectors, 
with an additional window to account for uncertainty in the measurement of the time of 
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arrival. The triggers must also be observed with the same template in both detectors. 
Triggers that survive the time and parameter coincidence test are referred to as candidate 


events. Details of the coincidence tests are presented in Section 3.5 


The quadrature sum of the reweighted SNR p in each detector is the pipeline’s 
detection statistic which ranks the likelihood that a trigger is due to a gravitational-wave 
signal. To assign a statistical signihcance to detection candidates, the pipeline measures 
the false-alarm rate of the search as a function of the detection-statistic value pc- Since 
it is not possible to isolate the detectors from gravitational waves, it is impossible to 
directly measure the detector noise in the absence of signals. This, together with the 
non-stationary and non-Gaussian nature of the noise, means that the false-alarm rate 
of the search must be empirically measured. This is done by applying a time shift to 
the triggers from one detector relative to another. The minimum time-shift offset is 
chosen to be larger than the time-coincidence window used to determine if signals are 
observed with consistent parameters in the network. Events in the time-shifted analysis 
therefore cannot be due to the coincidence of the pair of triggers produced by a real 
gravitational-wave signal. Many time shifts create a large background data set which 
are used to approximate the background noise and estimate the search’s false-alarm 
rate. 


Since different templates in the bank can respond to detector noise in different 
ways, the search background is not uniform across the template bank. To maintain 
the sensitivity of the search to signals over a wide range of masses under this non- 
uniform background distribution, the search sorts both candidate events and background 
events into different classes. The false-alarm rate of the search in each class is used to 
assign a p-value to the candidate events; a given candidate event is compared to the 
background events from the same class. To account for having searched multiple classes 
whose response is not completely independent [HS], the significance of candidate events 
is decreased by applying a trials factor equal to the number of bins, ribins to obtain a final 
p-value which describes the statistical significance of a candidate event. This procedure 
is described in Section 13.61 


3. Search Description 

In this section we describe the methods and algorithms used in the PyCBC search 
pipeline introduced in Section In particular, we focus on the parts of the PyCBC 
search that improve upon the ihope pipeline described in Ref. |18j . 

3.1. Template Bank Placement 

Methods for creating template banks of non-spinning and aligned-spin waveform hlters 
for a given parameter space have been extensively explored in the literature |l2l 0310H 
05l06li7l08l09l[5Ol|5ll[52l[53ll6l]. Here we use the method presented in [53]. The 
density of templates across the parameter space depends upon the noise power spectrum 





PyCBC search 


of the detectors. 

The PyCBC search generates a template bank covering the four dimensional space 
of mass and aligned spin of the two components. The pipeline uses a single template 
bank for all detectors and over the entire duration of the search. Using a common 
template bank for all templates allows us to require coincident events to be observed 


in the same template in different detectors, as discussed in Section |3.5[ which improves 
the sensitivity of the pipeline. 

To generate an appropriate bank, we must compute a noise PSD estimate averaged 
both over time and over detectors. We have explored several different methods to 
create a single noise PSD that is valid for the duration of the analysis period [HS] and 
hnd that the harmonic mean provides the best estimate for placing the template bank. 
Specihcally, we measure the noise PSD every 2048 seconds over the observation period, 
independently in each detector, using the median method of Ref. HT]. We then obtain 
Ng power spectra S'„ for each detector in the network. We hrst construct the harmonic 
mean PSD for a single detector, dehned by averaging each of the fk frequency bins 
independently according to 




harmonic 


ifk) = N, 



(7) 


We then use the same method to compute the harmonic mean of the resulting PSDs 
from each detector in the network. 

The PSD estimate, and hence the bank itself, only needs to be re-generated when 
there are signihcant changes in the detector’s noise PSD. This typically only happens 
when signihcant physical changes occur in the detectors. Using a single PSD estimate 
for an extended period of time allows the efficient use of template banks that include 
compact-object spin, as demonstrated in Ref. [2H]- In Section 4, we compare the 
sensitivity and the computational cost of using a single template bank constructed 
with an averaged PSD and that of using regenerated template banks with shorter PSD 
samples. 


3.2. Removal of Non-Gaussian Noise Prior to Filtering 

Transient noise in the detectors’ data streams can produce high SNR triggers, even 
when the transients do not resemble the templates. Data quality investigations and 
vetoes gD [33] remove many, but by no means all, of these loud transients which 
can affect the astrophysical sensitivity of the search. Although short-duration, loud 
transients (or glitches) that survive data-quality investigations are suppressed by the chi- 
squared test, they can reduce search sensitivity through two mechanisms: dead time due 
to the clustering algorithm and ringing of the matched hlter; both of these mechanisms 
are related to the impulse response of the matched hlter. In this section, we explain the 
origin of these features and describe the methods used by the pipeline to reduce their 
effect on the search’s sensitivity. 
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Figure 2. The effect of gating on a loud noise transient. Both lines show the detector 
strain, which has been rescaled by a factor of 10^^ prior to filtering. On this scale, 
the transient has a peak magnitude over 5,000. The blue line shows the data before 
applying the Tukey window, the red line shows the data after. 


The impulse response of the matched hlter is obtained by considering a delta- 
function glitch in the input data s{t) = 6{t — tg), where tg is the time of the transient. 
Although not all types of glitches are of this nature |66] , many loud noise transients are 
well approximated as s{t) = n{t) + 5{t — tg). For example, Figure]^ shows a typical loud 
transient glitch from LIGO’s sixth science run. As loud noise transients occur frequently 
in detector data, most glitches are not individually investigated. Although the cause 
of this glitch is unknown, it likely comes from an instability in the detector’s control 
systems. For such a glitch, the matched-hlter SNR given by Eq. ([^ will be dominated 
by 


P^{t) ~ llosit - tg) + Iln{t - tg), 


( 8 ) 


where 

IcosMt - tg) = 4Re hco.M f) 

'-'/low ) 

is the impulse response of each of the two phases of the matched hlter. This is equal to 
the convolution of the template, hcos,sm{t) with the inverse Fourier transform of inverse 
power spectral density, 1/S'„(/). To ensure that the impulse response of the hlter is 
of hnite duration, the inverse PSD is truncated to a duration of 16 seconds in the 
time domain before hltering HT]. The inverse PSD truncation serves to smear out very 
sharp features in the PSD: using a 16-second window provides sufficient resolution of 
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Figure 3. Response of the pipeline to a loud glitch. The matched-filter SNR of 
the triggers generated is shown as a function of time. Black circles show the triggers 
generated immediately after filtering without gating applied to the input data. Blue 
triangles show the triggers that remain after the re-weighted SNR is computed for each 
trigger and a threshold of p > 5.5 is applied. Although the significance of these triggers 
is suppressed by the re-weighted SNR, many triggers still remain which can increase 
the noise background of the search. The red crosses show the triggers produced by 
the search after gating is applied and the SNR is reweighted. The large majority of 
triggers caused by the glitch in the ±8 seconds around the transient are absent. 


these features, while minimizing the length of the impulse response. The length of the 
template can be on the order of minutes, leading to very long impulse response times 
for the matched hlter, although the template amplitude is quite sharply peaked near 
merger. 

Figure l^shows the result of hltering the glitch in Figure [^through a template bank 
and generating triggers from the matched-hlter SNR time series. The black circles show 
the triggers generated by the ringing of the hlter due to the glitch. The largest SNR 
values occur close to the time of the glitch and are due to the impulse response of the 
template h. The shoulders on either side of this are due to the impulse response of 
the inverse power spectrum 1/Sn{f)- This leads to the two effects described above: an 
excess of triggers around the time of the glitch which can increase the noise background 
of the search and a window of time containing multiple noise triggers that can make it 
difficult to distinguish signal from noise. In Figure we also plot the triggers that have 
a re-weighted SNR above the threshold p > 5.5. Although the use of the chi-squared 
veto to construct the reweighted SNR suppresses the signihcance of these triggers, it 
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does not completely remove them from the analysis and an increased trigger rate around 
the time of the glitch remains. 

The PyCBC pipeline identihes non-stationary transients and removes them from 
the data, a process called gating. The data around short-duration noise transients is set 
to zero prior to matched hltering. To zero the data, the input data s{t) is multiplied 
by a window function centered on the time of the peak. A Tukey window is applied to 
smoothly roll the data to zero and prevent discontinuities in the input data. The Tukey 
window used has a shape parameter of a = 0.5, where a is the fraction of the window 
inside the cosine tapered region. The effect of this gating on the strain data is shown 
in Figure where the input data is zeroed for 1 second around the time of the glitch. 
Figure shows the output of the search when run on the gated data. In addition to 
removing the loud triggers with SNR ~ 1000 at the time of the glitch, gating removed 
the additional triggers with SNR ~ 10 generated by the impulse response of the hlter 
before and after the glitch. This improves search sensitivity by reducing the amount 
of data corrupted by the glitch to only the windowed-out data, and reduces the overall 
noise background by removing noise triggers that could possibly form coincident events. 
Gating is typically only applied to short-duration transients length of order 1 s. Longer 
duration transients that are identihed by data quality investigations are removed prior 
to analysis by the pipeline. 

This process requires the identihcation of the time of noise transients that will be 
removed by gating. These times can be determined by either data quality investigations 
or by a separate search for excess power in the input strain data (typically with a high 
SNR threshold) [HU |68]; this method was used in Ref. [I]. Alternatively, short-duration 
glitches that are not flagged by data quality investigations can be identihed by the 
pipeline. To do this, the pipeline measures the power spectral density of the input 
time series. The data is Fourier transformed, whitened in the frequency domain, and 
then inverse Fourier transformed to create a whitened time-series. The magnitude of 
the whitened strain data is computed and times exceeding a pre-determined threshold 
are identihed. Peaks that lie above the threshold are identihed by the time-clustering 
algorithm of Ref. mi. The data around these times is then gated using the procedure 
described above. 

3.3. Matched Filtering 

A core task of the search pipeline is to correlate the detector data against a bank of 
template waveforms to construct the matched hlter signal-to-noise ratio. The matched 
hltering used in our search pipeline is based on the FindChirp algorithm developed 
for use in the Initial LIGO/Virgo searches for gravitational waves from compact 
binaries mi. In this section we describe two improvements to the FindGhirp algorithm; 
changes to the noise power spectral density estimation and a new thresholding and 
clustering algorithm that identihes maxima in the signal-to-noise ratio time series. 

The matched hlter in Eq. ([^ is typically written in terms of continuous quantities. 
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e.g. s{t),h{t), and Sn{f)- In practice, the pipeline works with discretely sampled 
quantities, e.g. Sj = s{tj), where Sj represents the value of s{t) at a particular time 
tj. Similarly, Sn{fk) represents the value of the noise power spectral density at the 
discrete frequency fk- The input strain data is a discretely sampled quantity with a 
hxed sampling interval, typically At = 1/4096 s. Fourier transforms are computed using 
the Fast Fourier Transform (FFT) algorithm in blocks of length Tb = 256 seconds. For 
this block length, the number of discretely sampled data points in the input time-series 
data Sj is N = Ts/At = 256 x 4096 = 2^°. The discrete Fourier transform of sj is given 
by 

N-l 

4 = ( 10 ) 

j=0 


where k = fk/ {NAt). This quantity has a frequency resolution given by A/ = l/(iVAt). 

To construct the integrand of the matched hlter, the discrete quantities hk, Sk, 
and Snifk) must all have the same length and frequency resolution. The most 
straightforward way to do this is to couple the computation of Sk and Sn{fk), by using the 
same length of data to compute both. An average PSD is then computed by averaging 
(typically) 15 overlapping blocks of data [69l [171 HH], leading to a large variance in the 
estimated PSD. 

The PyCBC pipeline decouples the computation of Sn{fk) and Sk, allowing many 
more averages to be used to compute Sn{fk)- As discussed in Section 3.2, the inverse 
PSD is truncated to 16 seconds to reduce the length of the hlter’s impulse response: 
thus, it is natural to estimate the PSD using 16-second time-domain blocks. An input 
hlter length of 2048 seconds allows for 127 such blocks, leading to signihcantly less 
variance in the PSD estimate. The resulting PSD is then linearly interpolated to 
match the resolution of the data. Finally, the interpolated Sn{fk) is inverted, inverse 
Fourier transformed, truncated to a hxed length of 16 seconds in the time domain im, 
and then Fourier transformed to the frequency domain for use in the matched hlter. 
Implementation of the matched hltering algoritm in PyCBC can be performed using 
either single-threaded or parallel FFT engines, such as FFTW 123 or the Intel Math 
Kernel Library. The choice of single- or multi-threaded hltering is made at runtime, so 
that the fastest implementation (one multi-threaded process, or several single-threaded 
processes) can be made depending on the architecture used. 

Once the matched-hlter SNR time series has been computed, the hnal step of the 
hltering is to generate triggers. These are maxima where the SNR time series exceeds a 
chosen threshold value. For either signals or noise transients, many sample points in the 
SNR time series can exceed the SNR threshold. Since a real signal will have a single, 
narrow peak in the SNR time series, the pipeline applies a time-clustering algorithm to 
keep local maxima of the SNR time series that exceed the threshold. 

The PyCBC pipeline divides the 256-second SNR time series into equal 1-second 
windows and then identihes the maximum of the time series within each window. 
As the maximization within each window is independent from other windows, the 
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clustering for each window can be parallelized for efficiency on multi-core processing 
units. If the desired clustering window is greater than 1 second, then the initial list of 
maxima—potentially one for each window within the segment—is further clustered: a 
candidate trigger in a window is kept only if it has a higher SNR than both the window 
before and after it, ensuring that it is a local maximum in the SNR time series. This 
clustering algorithm is more computationally efficient than the running maximization 
over template length described in Ref. HT]. Furthermore, it reduces the clustering 
window, thereby increasing the search sensitivity. 

3 . 4 . Signal Consistency Tests 

The chi-squared signal-consistency test introduced in Ref. im and developed in Ref. [I6] 
is a powerful way to distinguish between signals and noise in searches for gravitational 
waves from compact-object binaries. The PyCBC search allows the number of bins p in 
Eq. ([^ to be varied as a function of the parameters of the template [28]. This results 
in an improvement in search sensitivity, as demonstrated in Section 

The chi-squared statistic can be computationally intensive to calculate, however. 
For every FFT operation required to compute the SNR time series, the chi-squared 
statistic requires an additional p FFT operations, one for each frequency bin of Eq. ([^. 
The PyCBC pipeline instead calculates the chi-squared test only at the time samples 
corresponding to clustered triggers in the matched-filter SNR time series. To do this, 
the pipeline uses an optimized integral over frequency, rather than computing the FFT. 
If data quality is poor and the number of triggers in a FFT block is large, the FFT 
method becomes more efficient. We quantify this cross-over point below. 

To compare the computational cost of the two methods, we consider the calculation 
of the p matched filters for the naive chi-squared test implementation. Computation of 
the pf for each of the p bins requires an inverse complex FFT. For a data set containing 
N sample points the number of operations is p x 5A^log(iV). [|] If we instead calculate 
the chi-squared value only at the peaks, we must evaluate 



( 11 ) 


at the set of points [j] identified as triggers. We can re-write Eq. (0 as 



( 12 ) 


where denote the frequency boundaries of the p bins and are given by fc" 

jmm,ma,x^^ and Qk is the kernel of the matched filter, defined as [69] 


mm,max 


mm,max 


mm,max 


- Hf)Hf) 
S„(f) 


(13) 


f Since the FFT dominates the operations count, we have neglected lower-order terms that do not 
significantly contribute to the computational cost. 
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We can further reduce the computational cost of the chi-squared by noting that the 
exponential term requires the explicit calculation of at fcmax = points. 

This can be reduced to a single computation of the exponential term by pre-calculating 
g27rij/Af then iteratively multiplying by this constant to obtain the next term 


needed in the sum. To do this, we write Eq. (12) in the following form: 


V 


kT 


iJi = E E 




2TTij /N'^ 


k-1 


(14) 


2 = 1 




This reduces the computational cost of each term in the sum to two complex 
multiplications: one to multiply by the pre-computed constant and one for 


the multiplication by q. Computing the right hand side of Eq. (14) then requires the 


addition of two complex numbers for each term in the sum. The total computational 
cost to compute the chi-squared test for Np points is then 14fcmax x Np = 7N x Np. 

For small numbers of matched-hlter SNR threshold crossings, this new algorithm 
can be significantly less costly than calculating the chi-squared statistic using the EFT 
method. However, if the number of threshold crossings Np is large, then the EFT 
method will be more efficient due to the logA^ term. The crossover point can be 
estimated for p chi-squared bins as 


px5N\og(N) 5 , 

Uk -7P‘°sW. 


(15) 


although this equation is approximate because the computational cost of an EFT is 
highly influenced by its memory access pattern. For a typical LIGO search where 
N = 2^^, the new algorithm is more efficient when the number of points at which the 
statistic must be evaluated is Np < 100. For real LIGO data, the number of times that 
the statistic must be evaluated is found to be less than this threshold on average, and 
so this method significantly reduces the computational cost of the pipeline. However, 
there are still periods of time where the data quality is poor and the FFT method is 
more efficient. Gonsequently, the pipeline described here computes Np and uses either 


the single-trigger or FFT method, depending on the threshold determined by Eq. (15). 


For a typical analysis configuration that only identifies triggers once every second, this 
threshold is never exceeded, and the single-trigger method is faster even where the data 
quality is poor. 


3.5. Coincidence Test 

To reduce the rate of false signals, we require that a candidate event is seen with 
consistent parameters in all of the detectors in the network. The pipeline enforces 
this requirement by performing a coincidence test on triggers produced by the matched 
filter. This test is performed after triggers occuring during times affected by known 
instrumental or environmental artifacts have been discarded. The PyGBG pipeline 
requires consistency of arrival time and template parameters (masses and spins) between 
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different detectors. Here, we only consider a two-detector network, but the extension of 
this coincidence test to more than two detectors is straightforward. 

For the two-detector LHO-LLO network, signals must be seen in both detectors 
within a time difference of 15 ms: 10 ms maximum travel time between detectors 
and 5 ms padding to account for timing errors. Since the same waveform should be 
observed in all detectors, the pipeline requires consistency of the parameters of the 
best-fit template between detectors. Previous searches used a metric-based coincidence 
test 1721 that checked the consistency of the template parameters and arrival times of 
triggers between detectors, but did not require template parameters to be the same. 
The PyCBC pipeline requires an exact-match coincidence. The same template bank is 
used to filter data and produce triggers from each detector. In addition to enforcing 
coinsistent arrival times, the intrinsic parameters (masses and spins) of triggers in each 
of the detectors must be exactly the same. The exact-match coincidence test is useful 
in cases where there is no simple metric to compare gravitational waveforms, such as 
template waveforms for binaries with spinning neutron stars or black holes ITS] or for 
high-mass waveforms where a stochastic template placement algorithm is used [53l |28] . 
As demonstrated in Section]^ the use of exact match coincidence leads to a measurable 
improvement in the search sensitivity, even when a metric-based coincidence test is 
available. 

Triggers that survive the coincidence test are considered coincident events. These 
candidates are then ranked by the quadrature sum of the reweighted SNR in each 
detector; for a two-detector network, we have 



( 16 ) 


The choice of coincidence test has implications for the number of single-detector 
triggers that must be stored. In previous searches, triggers were clustered over the entire 
template bank [ZZIIIH] prior to generating coincident events. However, applying such 
clustering in conjunction with exact-match coincidence leads to unacceptable reductions 
in search sensitivity, as random noise causes the highest-SNR template to vary between 
detectors for quiet simulated signals. Thus, we keep all triggers from the bank before 
testing for coincidence. 

3.6. Candidate Event Significance 

The final step of the pipeline is to measure the false-alarm rate of the search as a function 
of the detection statistic, pc, and use this to assign a statistical significance to candidate 
events. The rate and distribution of false alarms in the search, i.e. coincident triggers 
due to noise, depends on the pipeline’s response to unpredictable, non-Gaussian and 
non-st at ionary detector noise and must be measured empirically. The PyCBC pipeline 
measures the false-alarm rate using time shifts. In this section, we describe the method in 
detail, again restricting attention to the two detector case. Generalisations to additional 
detectors are relatively straightforward. 
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Triggers from one detector are shifted in time with respect to triggers from the 
second detector, then the coincidence test is re-compnted to create a backgronnd data 
set that does not contain coincident gravitational-wave signals. Repeating this procednre 
many times on a snfficiently large dnration of inpnt data prodnces a large sample of false 
coincidences that are used to compute the false-alarm rate of the search as a function 
of the detection statistic. Under the assumptions that the times of transient noise 
artifacts present in the data streams are not correlated between different detectors and 
that gravitational-wave signals are sparse in the data set, this is a sufficiently accurate 
approximation of the background to support detection claims. 

Our signihcance calculation requires that candidate events are statistically 
independent; however, both noise transients and signals can generate many triggers 
across the template bank, potentially yielding several correlated coincident events within 
a short time (< 1 second). To generate independent candidates, the pipeline performs a 
hnal stage of clustering: if more than one coincident event occurs within a time window 
of hxed duration, typically 10 seconds, only the event with the highest detection statistic 
value An identical clustering operation is also performed on the events in each time- 
shifted analysis. 

Each candidate event is assigned a p-value that measures its signihcance. For 
a candidate event with detection statistic Pc, its p-value pb is the probability that the 
there are one or more coincident noise events (false alarms) that have a detection statistic 
value greater than or equal to pc- We calculate p-values under the null hypothesis that 
all triggers seen are due to noise. While this hypothesis might be unrealistic in the 
presence of loud signals, note that we can never determine with complete certainty from 
gravitational-wave data alone that any given trigger is due to signal rather than noise. 
In order to claim detection it is necessary to show that the statistic values of events 
actually obtained in the search are highly improbable under such a null hypothesis: i.e. 
that the p-value is very small. As demonstrated in m, including all search triggers in 
the background calculation results in a self-consistent signihcance, i.e. small p-values are 
assigned to random noise candidate events with the expected frequency; the procedure 
also does not adversely affect the efficiency of the search compared to other methods 
which aim to remove likely signal triggers before calculating the signihcance of search 
events. 

Using the distribution of coincident events from the time shifts, we can measure 
how many noise background events rib are louder than a given candidate event. We thus 
determine the function nb{pc) which gives the number of background events having a 
higher detection-statistic value than p^. The probability that one or more noise events 
as loud as a candidate event with detection-statistic value p* occurs in the search, given 
the duration of observing time T and the amount of background time constructed from 
the time shifts Tb, is 


p(> 1 above pl\T, Tb)o = 1 — exp 


-T{l + nb{pl)) 
Tb 


A detailed derivation of this result is given in Appendix A 


( 17 ) 
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To account for the search background noise varying across the target signal space, 
candidate and background events are divided into three search classes based on template 
length. The significance of candidate events is measured against the background from 
the same class. This binning method prevents loud background events in one class 
from suppressing candidates in another class. For each candidate event we compute 
the p-value inside the class A, p^{pc) = 1 — exp [—T(1 + n^{pc))/T}^ , where n^{pc) 
is the number of noise background events in the class A above the candidate event’s 
re-weighted SNR value pc- To account for having searched multiple classes which could 
all give rise to false alarms, the significance is decreased by a trials factor equal to the 
number of classes nc [M]; the final significance of a candidate found in class A is given 
by 


PbiPc] A) 


1 — exp 


-ncTjl + n^jpc)) 
Tb 


(18) 


Equivalently, we may consider the quantity nc(l -l- n^{pc))/Tb as a rate of false positive 
events more signihcant than such a candidate, when summed over the search as a whole. 

The observation time used in the search Tobs is determined by two considerations. 
The maximum length of the observation time is set by the requirement that the data 
quality of the detectors should not change significantly during the search. If the time- 
shifts mix data from periods of substantially different data quality (i.e. higher- or lower- 
than-average trigger rates), the noise background may not be correctly estimated. The 
minimum length of the observation time is set by the false-alarm rate necessary to 
claim the detection of interesting signals. The smallest false-alarm rate that the search 
can measure scales is achieved by performing every possible time shift of the two data 
streams. While the length of data analyzed in any given time shift will vary, the total 
amount of background time analyzed will equal Tf, = T^/5, where 6 is the time-shift 
interval. [|] Thus, the minimum false-alarm rate scales as 6/T^^^. In a two-detector 
search using time shifts of 0.1 seconds, approximately fifteen days of coincident data 
are sufficient to measure false-alarm rates of 1 in 2 x 10® years, corresponding to a 
significance of 5 a. This is the length of data used in the observation of GW150914 [1]. 

Two computational optimizations are applied when calculating the false-alarm rate 
of the search. The implementation of the time-shift method used does not require 
explicitly checking that triggers are coincident for each time shift successively. Instead, 
the pipeline takes the triggers from a given template and calculates the offset in their 
end times from a multiple of the time-shift interval. It is then possible to quickly find the 
pairs of triggers that are within the coincidence window. Furthermore, the number of 
background triggers is strongly dependent on the detection-statistic value. For Gaussian 
noise, the number of triggers will fall exponentially as a function of pc- At low detection- 
statistic values, it is not necessary to store every noise event to accurately measure the 
search false-alarm rate as a function of detection statistic. Instead, the pipeline is given 


§ In the case where gaps between consecutive portions of data are a multiple of the time-shift interval 
this is exact. This is generally the case in our analyses. If gaps are not multiples of the time-shift 
interval, the total analyzed time will be approximately Tf, = T^/d, but there will be small corrections. 
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Time (days since Jun 26 2010 00:00:00 UTC) 


Figure 4. Sensitivity of the gravitational-wave detectors for the last part of the sixth 
science run for LIGO (S6D) and the third VIRGO science run (VSR3). The plot shows 
the volume-weighted average distance at which a 1.4, 1.4 BNS would be observed with 
an signal-to-noise ratio of 8 for each detector. The two rectangles indicate time intervals 
used for this study, corresponding to a week of data from July 2010 and a week of data 
from August 2010. 


a threshold on pc and a decimation factor d. Below this threshold the pipeline stores 
one noise event for every d events, storing the decimation factor so that the false-alarm 
rate can be correctly reconstructed. This saves disk space and makes computation of 
the false-alarm rate at given values of the detection statistic faster. 


4. Comparison to Initial LIGO Pipeline 


In this section, we compare the PyCBC search pipeline to the ihope pipeline used in 
the previous LIGO and Virgo searches for compact-object binaries [18]. We focus on 
tuning and testing the pipeline to improve the sensitivity to binary neutron star systems, 
although we note that use of the pipeline is not restricted to these sources. This pipeline 
has been used to search for binary black holes m [28113I, binary neutron stars and 
neutron-star black-hole binaries na in Advanced LIGO data. Section |4.1| describes 


the method that we use to measure the sensitivity of this pipeline and Section [4.2 
compares this sensitivity to that of the ihope pipeline. The comparison is performed 
using two weeks of data from the sixth LIGO science run [271133] shown in Figure 
We demonstrate that for binary neutron stars, the PyCBC search can achieve a volume 
sensitivity improvement of up to 30% over the ihope pipeline without reducing the 
sensitivity to higher mass systems. 
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As the metric of search sensitivity, we measure the pipeline’s sensitive volume. This is 
proportional to the fraction of sources that the pipeline can detect per unit time at a 
given false alarm rate if, given by 

V{J^) = f e(J^; X, A)0(x, A)dxdA. (19) 


Here, A are the physical parameters of a signal, x is a spatial co-ordinate, 0(x, A) is the 
distribution of signals in the universe, and e is the efficiency of the pipeline at detecting 
signals with parameters A in volume x and with false alarm rate iF. 

We find the sensitive volume by adding a large number, Nj, of simulated signals 
to the data and measuring the pipeline’s ability to identify them. If the simulated 
signals are drawn from the same distribution as the astrophysical distribution 0, then 
the sensitive volume is 

1 

^ i=i 

where W, is the false alarm rate associated to simulation i and 0(W|Wj) is a step function 
with 0(W|Wi) = 1 if Wi < J- and zero otherwise. The false alarm rate for a simulation 
is calculated using the most significant event within a 1 second window of the arrival 
time of the simulated event. If no event is found within this time, the detection statistic 
value is set to zero. 

For this comparison, we assume an astrophysical distribution 0 in which signals 
are uniformly distributed in volume, and isotropic in sky location and orientation, [jj] 
Distributing simulations uniform in volume leads to a majority of signals being injected 
at large distances, where almost all will not be detected due to less than optimal sky 
position or orientation. For a fixed number of injections, this produces large errors 
in the estimate of V. Instead, we generate signals distributed uniformly in distance 
between rmin and rmax chosen such that all signals with distance r < rmin will be found, 
even at small (< 10“^/yr) false alarm rate thresholds, and all signals with r > r^ax 
will be missed, even at large (> 10^/yr) W. Gravitational waves from more massive 
systems have larger amplitudes and can be detected at greater distances than less 
massive systems. The amplitude of the emitted signal scales, at leading order, with 
the chirp mass of the binary Ai = /{mi + 7712 )^^®. If a binary with chirp 

mass M.Q can be detected at a distance of tq, then a binary with chirp mass M. can be 
detected at a distance given approximately by 

D = rii{MJMof'^. (21) 


Consequently, we use mass-dependent bounds ryam,i and rmax,i, scaled as in Eq. (21), 
for a simulated signal with chirp mass ATj. The sky locations and orientations of the 


II At small distances, it is necessary to account for discreteness of galaxies while, at larger distances, 
cosmological effects become important. 
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simulate population are distributed isotropically. Since the distribution of simulated 
signals differs from the assumed astrophysical distribution 0, we must re-weight the 
contribution each simulation makes to the sensitive volume as [53] : 


Nl r 




i=l 




= ■ 


where Ar* = The error on this estimate is 


( 22 ) 


sv = (23) 

To quantify the sensitive volume of each pipeline for this study, we select a 
compact-object binary population distributed uniformly over component masses with 
1 A < 7 and mi-|-m 2 < 14 Mq. The parameters A are therefore the component 

masses of the binary mi, m 2 and the orientation and polarization angle of the binary 
with respect to the detectors. For data from LIGO’s sixth science run, we hnd that 
?"min = 0.5 Mpc and r^ax = 30Mpc are suitable choices for a binary with component 
masses mi = m 2 = 1.4 M©. 


4 . 2 . Relative Search Sensitivity and Computational Cost 

Both the PyCBC and ihope [H] pipelines were used to analyze two weeks of data 
from LIGO’s sixth science runs and the sensitive volume of each search, as well as the 
computational cost, was compared. Both searches use a template bank designed to 
search for compact-object binaries with component masses between 1 and 24M©, total 
mass mi -|- m 2 < 25 Mq and zero spins, as shown in Fig. In the literature it is 
generally assumed that the noise power spectral density is invariant, which is not the 
case for real data. In searches of initial LIGO and Virgo data, the time-dependence of the 
detector PSDs was addressed by recomputing the template bank on regular intervals of 
around an hour [18]. Since PyCBC’s templates include both mass and spin parameters, 
regenerating the bank so often is not computationally feasible; instead, the use of a 
single, pre-generated bank is required, as described in Section |3.1[ 

The PyCBC analysis was conhgured as described in detail in Sections]^ andThe 
chi-square test uses p = 100 bins for templates with a chirp mass Ai < 1.74 Mq and 
p = 16 for higher masses. To estimate the noise background of the search, an 0.2 s 
time-shift interval is used, and events are divided into two classes: one class contains 
triggers with M. < 1.74 Mq, a second triggers with 1.74 < A1/Mq < 6.1, and templates 
outside these bins are ignored in the search. 

The ihope analysis was conhgured as for the S6-VSR2/3 search for low-mass 
binaries [2Z|. To reduce computational cost, the ihope pipeline has two stages of matched 
hltering, where the computationally costly chi-square test is performed at the second 
stage on only those triggers found in coincidence (and time-shift coincidence). When 
large numbers of time shifts are performed, the majority of triggers appear in coincidence 
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Figure 5. Mass-ranges for software injection, shown in the mi — m 2 mass-plane. The 
template bank used to search for these injections is indicated by hatched regions and 
the injection set by the red shaded region. The black dashed lines show chirp masses 
of 3.48 Mq and 6.1 Mq, the boundaries between the mass bins used. Triggers from 
templates with chirp masses larger than 6.1 Mq are discarded in post-processing. 


in at least one time shift and there is little beneht to the two-stage analysis. Thus, to 
investigate sensitivity at low false alarm rates, we also run the pipeline with a single 
stage of matched hltering where the chi-squared test is performed on every trigger. 

Figures and compare the sensitivity of the two pipelines, using a week of data 
from July 2010 and a second week from August 2010. The data from July had larger 
fluctuations in the overall sensitivity, but relatively few noise transients while the August 
data had a consistent range but more glitches [3S]- Using these two different weeks allows 
us to test the pipelines under various conditions. Figure [^restricts to low mass binary 
coalescence signals (AJ < S.JSMq) while Figure displays higher masses. The hgures 
show that the sensitive volume of the new pipeline is greater than or equal to that of 
the ihope pipeline for both subsets of signals in both weeks of data. 

At a false-alarm rate of 1 in 100 years, the PyCBC pipeline affords a sensitivity 
improvement for low-mass systems of 20% for the hrst week of data and 30% for the 
second week of data. The matched-hlter SNR of simulated signals does not change 
signihcantly between the old and new pipelines. Furthermore, the chi-squared value, 
and hence re-weighted SNR, is also similar for both pipelines. The improvement in 
sensitivity of the PyCBC analysis is due to the reduction in the noise background of the 
search. Approximately 10% of the improvement in sensitive volume comes from using a 
hxed template bank with the exact-match coincidence test. The remaining improvement 
is due to reduction in the noise background because of the increased number of chi- 
squared bins used in the low-mass region of the search. Both of these improvements 
allow us to decrease the number of high signihcance events caused by noise transients. 





Volume (Mpc^) Volume (Mpc^) 


PyCBC search 


22 



55000 


50000 


45000 


40000 


35000 


S6 Two Stage 
Advanced Search 
■ “ S6 Single Stage 


False Alarm Rate (yr ') 



Figure 6. The sensitive volume for low-mass binary systems, chirp mass At < 3.48M0, 
of the PyCBC and ihope searches as a function of the false alarm rate threshold. The 
left plot shows results from a week of data from July 2010 while the right plot uses 
data from August 2010. The PyCBC search (red, solid line) is more sensitive over a 
broad range of false alarm rates than the ihope single stage (blue, dashed line) and 
ihope two-stage (pink dotted line) searches. The error bars for the ihope analyses are 
larger than for PyCBC as fewer injections were performed. Since it uses only 100 time 
shifts, the two-stage ihope analysis can only determine false alarm rates as low as six 
per year in one week of data. 



Figure 7. The sensitive volume for high-mass binary systems, chirp mass A1 > 
3.48M0, of the PyCBC and ihope searches as a function of the false alarm rate 
threshold. The left plot shows results from a week of data from July 2010 while 
the right plot uses data from August 2010. The PyCBC search (red, solid line) is more 
sensitive over a broad range of false alarm rates than the ihope single stage (blue, 
dashed line) and ihope two-stage (pink dotted line) searches. The error bars for the 
ihope analyses are larger than for PyCBC as fewer injections were performed. Since it 
uses only 100 time shifts, the two-stage ihope analysis can only determine false alarm 
rates as low as six per year in one week of data. 
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The improvement is most apparent in the August data as the rate of non-Gaussian 
noise transients in that week is larger than in the July data. Consequently, the new 
pipeline’s reduction of the noise background has a larger effect on the sensitivity for 
the August 2010 data. For high-mass systems, the PyCBC pipeline is marginally more 
sensitive than the ihope pipeline, with a 5-10% improvement at a false-alarm rate of 1 
in 100 years. The improvement is less pronounced as the number of chi-squared bins 
was unchanged at higher masses. A discussion of tuning the pipeline to binary black 
hole systems can be found in Ref. [28] . 


Job Type 

Two-Stage ihope 

Single-Stage ihope 

PyCBC 

Template Bank Generation 

13.3 

13.3 

4.7 

Matched filtering and 

515.4 

1431 

515.5 

Second Template Bank 

0.1 

- 

- 

Coincidence Test 

0.3 

8.3 

9.9 

Total 

529.1 

1453 

530.0 


Table 1. The computational costs of different parts of the single-stage and two-stage 
ihope search pipelines, and the new PyCBC pipeline. The costs are given in CPU days. 

Table [T] shows the computational cost of the pipelines. The cost of the PyCBC 
search is comparable to that of the two-stage ihope search, and about one third of the 
cost of the single-stage ihope search. However, the two-stage search is unable to compute 
false-alarm rates low enough to support a detection claim and is therefore unsuitable 
for analysis of advanced detector data. 

5. Conclusions 

This paper presents a new pipeline to search for gravitational waves from compact 
binaries developed in the PyCBC framework. The pipeline includes several new 
developments, including: (i) gating to mitigate the effect of loud, non-Gaussian noise 
transients in the input data; (ii) the ability to use a single, fixed template bank for the 
entire analysis created using the harmonic mean of a large number of noise spectral 
density estimates; (iii) decoupling of the matched filtering and noise power spectral 
density estimation to improve the performance of the matched filter; (iv) a simpler, more 
efficient algorithm for identifying and clustering peaks in the matched-hlter SNR time- 
series; (v) a computationally-efhcient implementation of the chi-squared veto; (vi) the 
ability to set the number of bins in the chi-squared veto as a function of the parameters 
of the template; (vii) a stricter coincidence test that requires that a trigger is found 
with exactly the same template in both detectors; and (viii) improvements to the time- 
shift method used to measure the false-alarm rate of the search and assign a p-value 
to candidate events. The new pipeline provides an improvement in sensitivity over the 
ihope pipeline na used in searches of initial LIGO and Virgo data 123. with up to 30% 
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increase in sensitivity to binary neutron star sources. The PyCBC pipeline described 
here was used in the LIGO Scientihc Collaboration and the Virgo Collaboration’s search 
for binary black holes in the hrst science run of Advanced LIGO [2HI El |T2], including 
the identihcation of the binary black hole signals GW150914 [1] and GW151226 |2]. 

The PyGBG framework provides a modular set of tools and algorithms that can be 
used to create flexible, sensitive gravitational wave search pipelines. For comparison 
presented in Section we deliberately conhgured the pipeline to mimic the ihope 
analysis as closely as possible. However, the PyGBG infrastructure can be conhgured 
in many different ways. For example, the mass range can be signihcantly extended 
[2H1 E], and the choice of template waveform can be varied across the mass space: e.g. 
Refs. [U |28] use post-Newtonian templates for system with a total mass less that 4 Mq 
and templates from the Effective One Body family for higher-mass systems PlEniES]. 
In principle, it is also possible to use the methods described here to implement matched 
hltering to detect gravitational waves from modeled transient systems other than 
compact-object binaries. 

Further improvements to the astrophysical sensitivity of the PyGBG search pipeline 
are ongoing. These include the extension of the pipeline to multiple detectors, upgrades 
to the detection statistic to include signal amplitude and phase consistency between 
detectors and the incorporation of methods to better measure the variation of the 
noise event distribution across the template bank. Ongoing advancements to the 
computational efficiency include investigating changes to the FFT block size and input 
data sample rate and the implementation of hierarchical search methods in the matched 
hlter. The FFT engine used in this pipeline can be efficiently used on graphics processing 
units (GPUs) and work is ongoing to develop a fast, efficient GPU implementation of 
this pipeline. In the longer term, the pipeline will be extended to search for binaries that 
exhibit spin-induced precession [70], include the effects of higher-order multipoles [02], 
perform a coherent search using a global network of detectors EH and target the search 
to look for gravitational waves emitted at the time of short gamma-ray bursts or other 
astronomical transients [701170]. Results from these studies will be reported in future 
publications. 
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Appendix A. Calculation of Event Significance 

In this Appendix we give the derivation of Equation ( |I^ which gives the probability 
that one or more events with higher detection statistic value than a given threshold p* 
would occur due to random coincidences of noise triggers over the search duration^] 
For a randomly chosen candidate event which is assumed to be due to noise, the 
probability that there are n* or fewer background events with a higher value of p* is 
given by 



where N}, is the total number of background events and Ne is the number of candidate 


events under consideration. This result follows by comparing a random (noise) candidate 
event’s pc value to an ordered list of background events from the time shifts. The 
candidate value can lie above all background events, below all background events or he 
in between two background events. There are Nt + 1 places where a given candidate 
event can lie when ranked against the list of background events. If the candidate event 
is due to noise then its detection statistic value is drawn from the same distribution as 
the time-shifted background events, therefore it is equally likely to occupy any one of 
these Ai"fe -|- 1 positions. 

Now, since higher pc values correspond to smaller rib values, the probability that 
one random coincident noise event lies above a threshold p* is 



(A.2) 


The same count of louder background events rib might be obtained for a slightly smaller 
value of Pc given that background event values are not inhnitely hnely spaced, thus 
the condition pc > p* is more restrictive than rib < and so strictly this equation 

should have a < sign; however we will neglect this subtlety in what follows. We wish 
to hnd the probability that one or more events out of the Ne coincident events in the 
search is a noise event (false alarm) at or above the threshold p*. This is given by the 
complement of the probability that all events do not lie above this threshold, which for 
one event is given by 1 — [1 -|- nb{p*e)] / [1 + Nb]. The probability that this is true for all 
candidate events follows by multiplying the individual probabilities for each of the Ne 
events: 



(A.3) 


^ A less detailed derivation was presented in [73]. 
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This step requires that candidate events are independent, which is achieved, to good 
approximation, by the hnal event clustering described in Section |3.6[ Then the 


probability that at least one out of clustered candidate events is louder than a 
threshold p* (if all candidate events are due to noise) is 


p(> 1 above p*|iVe, N^)q = 1 


1 + »-fe(Pc 

1 + iV, 


Ne 


(A.4) 


In what follows we will approximate 1 + ~ A'";,, which is a negligible correction given 

that Nh is very large. 

Since we do not know and Ne before performing the analysis, we should treat 
these event counts as stochastic variables. However, since there is such a large number 
of background events, the statistical uncertainty on Ni, is negligibly small. We expect 
the rate of background noise events for the duration of background time analyzed to 
be equal to the rate of coincident events over the duration of foreground time analyzed. 


(A.5) 


T: 

W N, 

_ _ 

T ~ Th¬ 
under the assumption that the candidate events are all due to noise, we can model N^. 
as a Poisson distribution with a mean of {Ne)o = {T/Tb)Nb. 

We wish to know the signihcance of obtaining a candidate with a given p* 
without restricting to any specihc number of coincident events, therefore we marginalize 
Eq. (A.4) over W to obtain 


p(> 1 above p*|W)o = 1 above W)oP(A^e|A^fe)- (A.6) 

Ve 


Given that coincident noise events are approximated by a Poisson process that we 
measure using the time-shifted background events, we can End the unknown probability 
p{Ne\Nh). Letting the Poisson rate of coincident noise events be p = {Nf,T)/Th, then 
the probability of obtaining W events is 


p{Ne\Ni,) = p{N^\p) = p 


jv, exp(-p) 

NJ ■ 


(A.7) 


Substituting into Eq. (A.6) we obtain 


p(> 1 above p*|W)o = < 1 - 


Ve 


1 - 


1 + Mp*c 

N, 






jv,e xp(-p) 
W! ' 


(A.8) 


Since the sum of a Poisson distribution over all possible event counts Xlve jg 

unity, this simplihes to 


p{> 1 above p*|iVb)o = 1 - 






1 - 


1 + nb{pl) 

N, 


Ve 


exp(-p) 


W! 


(A.9) 
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Next, we multiply inside the summation by exp[—p(l —exp[/i(l —(equivalent 
to multiplication by unity), to rewrite Eq. (A.9|) as: 


p(> 1 above p*|iVb)o = 1 - 






1 + 

Nh 




exp 


X 






exp 


-fi{l + nb{pl)) 

Nb 


(A.IO) 


Setting fi equal to /i[l — (1 + Ub) /A"fc], we again identify a sum of the Poisson probability 
over all possible counts: 


p(> 1 above pl\Nb)o = 1 


.We exp (-A) 

—^ 71 —exp 


N, 


iv. 


-p{l + nb{pl)) 

Nb 


(A.ll) 


All but the last term in the sum total to one and we can re-write this using the Poisson 
rate p = T{Nb/Tb), giving 


p(> 1 above p*\T, Tb)o = 1 — exp 
which reproduces Eq. © in the main text. 


-Tjl + nbjp*)) 

Tb 


(A.12) 
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